System and method for processing images

ABSTRACT

A method of processing a plurality of time separated images comprises selecting a plurality of imaging units in each image; measuring a temporal difference in each imaging unit; and selecting temporal differences above a threshold limit.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a continuation of U.S. patent application Ser. No.13/122,379 to Lee filed on Apr. 1, 2011 and entitled “SYSTEM AND METHODFOR PROCESSING IMAGES”, which is a national stage of PCT/CA2009/001397filed on Oct. 2, 2009, entitled “SYSTEM AND METHOD FOR PROCESSINGIMAGES”, which claims the benefit of U.S. Provisional Application No.61/102,206 to Lee filed on Oct. 2, 2008 and entitled “SYSTEM AND METHODFOR PROCESSING IMAGES”, the content of which is incorporated herein byreference.

FIELD OF THE INVENTION

The present invention relates generally to image processing and, moreparticularly, to a system and method for processing image data obtainedby using a dynamic imaging technique

BACKGROUND OF THE INVENTION

Medical imaging encompasses techniques and processes used to createimages of the human body for clinical purposes, including medicalprocedures for diagnosing or monitoring disease. Medical imagingtechnology has grown to encompass many image recording techniquesincluding electron microscopy, fluoroscopy, magnetic resonance imaging(MRI), nuclear medicine, photoacoustic imaging, positron emissiontomography (PET), projection radiography, thermography, computedtomography (CT), and ultrasound.

Medical imaging can incorporate the use of compounds referred to ascontrast agents or contrast materials to improve the visibility ofinternal bodily structures in an image.

Medical imaging was originally restricted to acquisition of singularstatic images to capture the anatomy of an organ/region of the body.Currently, the use of more sophisticated imaging techniques allowsdynamic studies to be made that provide a temporal sequence of imageswhich can characterize physiological or pathophysiological information.

Dynamic medical imaging involves an acquisition process that takes many“snapshots” of the organ/region/body of interest over time in order tocapture a time-varying behaviour, for example, distribution of acontrast agent and hence capture of a specific biological state(disease, condition, physiological phenomenon, etc.). As the speed anddigital nature of medical imaging evolves, this acquisition data canhave tremendous temporal resolution and can result in large quantitiesof data.

Medical imaging technologies have been used widely to improve diagnosisand care for such conditions as cancer, heart disease, brain disorders,and cardiovascular conditions. Most estimates conclude that millions oflives have been saved or dramatically improved as a result of thesemedical imaging technologies. However, the risk of radiation exposurefrom such medical imaging technologies for patients must be considered.

In this regard, the increasing use of CT in medical diagnosis hashighlighted concern about the increased cancer risk to exposed patientsbecause of the larger radiation doses delivered in CT than the morecommon, conventional x-ray imaging procedures. This is particularly truewith the two-phase CT Perfusion protocol. For illustration purposes,using the dose-length product for a protocol provided by a commerciallyavailable CT scanner (GE Healthcare), the effective dose of a CT Strokeseries consisting of: 1) a non-enhanced CT scan to rule out hemorrhage;2) a CT angiography to localize the occlusion causing the stroke; and 3)a two-phase CT Perfusion protocol to define the ischemic region withblood flow and blood volume and predict hemorrhagic transformation (HT)with blood-brain barrier permeability surface area product (BBB-PS); is10 mSv, of which 4.9 mSv is contributed by the two-phase CT Perfusionprotocol. The CT Stroke series is predicted to induce 80 and 131additional cancers for every exposed 100,000 male and female acuteischemic stroke (AIS) patients, respectively, with the two-phase CTPerfusion protocol alone predicted to cause half of the additionalcancers (Health Risks from Exposure to Low Levels of Ionizing Radiation:BEIR VII. The National Academies Press, Washington D.C., 2006).

Unless the effective dose of the two-phase CT Perfusion protocol isreduced, the benefits of medical imaging will be undermined by theconcern over cancer induction. Furthermore, the concern over the risksof radiation exposure extends to other medical imaging techniques,particularly for pediatric patients or those patients undergoingrepeated exposure.

The use of statistical filtering techniques to increase the signal tonoise ratio in low dose radiation imaging has been described, forexample in U.S. Pat. No. 7,187,794 issued Mar. 6, 2007. However, atpresent the application of statistical filtering is limited toprojection data prior to the construction of images. Severaldisadvantages are associated with statistical filtering of projectiondata. For example, the computational burden when working with projectiondata is high because each image from a dynamic sequence is typicallyreconstructed from ˜1,000 projections. In addition, the filteringprocess has to be ‘hardwired’ into the image reconstruction pipeline ofthe scanner. As such, statistical filters of projection data lackflexibility as they are typically specific to a scanner platform.

Whether perceived or proven, concerns over the risk of radiationexposure will have to be addressed to assure patients of the long termsafety of medical imaging techniques. Accordingly, there is a need formedical imaging techniques that allow for a reduction in radiation dose.

It is therefore an object of the present invention to provide a novelsystem and method for processing images.

SUMMARY OF THE INVENTION

In an aspect, there is provided a method of processing a plurality oftime separated images comprising: selecting a plurality of imaging unitsin each image; determining a temporal difference for each imaging unit;and selecting temporal differences above a threshold limit.

In another aspect, there is provided a system for processing a pluralityof time separated images comprising: an interface for receiving aplurality of time separated images; and a processor for selecting aplurality of imaging units in each image, determining a temporaldifference for each imaging unit, and selecting temporal differencesabove a threshold limit.

In yet another aspect, there is provided a computer readable mediumembodying a computer program for processing a plurality of timeseparated images, the computer program comprising: computer program codefor selecting a plurality of imaging units in each image; computerprogram code for determining a temporal difference for each imagingunit; and computer program code for selecting temporal differences abovea threshold limit.

BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments will now be described, by way of example only, withreference to accompanying drawings in which:

FIG. 1 is a flow chart of an image processing method;

FIG. 2 schematically illustrates the positioning of imaging units in aplurality of time separated images;

FIGS. 3A to 3E show brain blood flow maps from CT images processed witha statistical filtering technique;

FIG. 4 graphically shows blood flow Figure of Merit calculated for theimages shown in FIGS. 3A to 3D using the regions outlined in FIG. 3E;and

FIG. 5 is a diagram of a system for implementing the image processingmethod shown in FIG. 1.

DETAILED DESCRIPTION

Dynamic medical imaging involves an acquisition process that takes many“snapshots” of an organ/region/body of interest (i.e. a target region)over time in order to capture a time-varying behaviour, for exampleuptake and/or wash out of a contrast agent. This acquisition processresults in the production of a plurality of time separated images. Themethod and system described herein involves processing of such timeseparated images.

Turning now to FIG. 1, the steps performed during processing of timeseparated images according to the method are illustrated. The pluralityof time separated images may be obtained from dynamic medical imaging ofa patient with or without the use of an injected contrast agent. Acontrast agent can selectively increase the contrast of the targetregion in an image, for example on the basis of the target region'sstructure or physiological state. For example, one may inject into apatient a compound which has a biophysical, molecular, genetic orcellular affinity for a particular organ, disease, state orphysiological process. Such contrast agents are selected to have aproperty that provides enhanced information to a given imaging techniqueby altering imaging conditions, for example by altering image contrast,to reflect the behaviour of the compound in the body. This may beachieved via increased X-ray attenuation at a localized site (e.g. forCT/X-ray), altered paramagnetic properties (e.g. for MRI) or the use ofa radioisotope for nuclear medicine/PET. Contrast agents for manyimaging techniques are well-known. In some cases contrast enhancementdoes not rely on an injected contrast agent, for example the use of“black blood” or “white blood” in magnetic resonance imaging (MRI) wherespecific pulse sequences are used to change the magnetic saturation ofthe blood and thus its appearance in the image, or tagged MRI sequenceswhich alter the magnetic behaviour of a particular tissue or fluid. Aninjected contrast agent, or a tissue or fluid with altered behaviour,may all be regarded as “imaging agents”.

The plurality of time separated images are not limited to any particularimaging technique and include, for example, dynamic medical imagingusing magnetic resonance imaging (MRI), computed tomography (CT),nuclear medicine (NM) or positron emission tomography (PET). Theplurality of time separated images are also not limited to particularimage types. For example, grayscale or colour images may be processed.

During the method, within each time separated image a plurality ofimaging units are selected (step 110). An “imaging unit” may be anydesired unit for separating an image into equivalent portions including,for example, a pixel, a plurality of pixels, a fraction of a pixel, avoxel, a plurality of voxels, a fraction of a voxel etc.

The imaging unit data is represented by a value, for example a digitalvalue, and can thus be quantified and measured. Typically, imageintensity is measured for corresponding imaging units in the timeseparated images. For dynamic imaging using contrast agents a temporaldifference can be determined with respect to contrast concentration. Thetemporal difference of contrast concentration can be represented as aplot of contrast enhancement as a function of time.

Two alternatives for determining temporal differences of contrastconcentration and selecting temporal differences are shown in FIG. 1. Inone alternative shown in FIG. 1, once the plurality of imaging units hasbeen selected from each time separated image at step 110, a temporaldifference is determined for each imaging unit (step 115). Thedetermined temporal difference or representations thereof are thenranked according to degree of observed temporal difference (step 120).Temporal differences above a threshold limit are then selected (step125). A processed image is then constructed on the basis of the selectedtemporal differences (step 150).

In the other alternative shown in FIG. 1, the temporal differences arenot ranked. Rather, the temporal data for imaging units are analyzed toselect temporal differences above a threshold limit using statisticaltechniques, and more particularly blind signal separation statisticaltechniques. For example, with respect to dynamic CT imaging using acontrast agent, once the plurality of imaging units have been selectedfrom each time separated image at step 110, a covariance matrix of theco-variations of all pairs of contrast enhancement plots in relation toa mean curve is established (step 130) and analyzed with a blind signalseparation statistical technique (step 135), such as with a principalcomponents analysis being applied to the covariance matrix.Eigenvector/eigenvalue pairs are calculated on the basis of thestatistical analysis (step 140). Resulting eigenvector/eigenvalue pairsabove a threshold limit are then selected on the basis of theireigenvalues (step 145). A processed image is then constructed using theselected eigenvectors, with different weightings of the selectedeigenvectors being used to construct individual imaging units (step160).

The plurality of time separated images may optionally be registered tocorrelate corresponding locations in the time separated images to eachother. Registration of the plurality of time separated images is onlyneeded to correct for gross misregistration. Otherwise, slight movementbetween the time separated images may be tolerated and may even beremoved by the processing method described herein.

As another optional step, the plurality of time separated images may bepreprocessed so that the time separated images are represented withquantifiable values or codes. For example, the time separated images maybe represented by digital values using known digitization techniques. Ifthe imaging data of the plurality of time separated images is providedin a digital form, then a digitization step is not necessary.

The determination of temporal differences of contrast concentration willnow be further described with reference to FIG. 2. FIG. 2 shows a seriesof six chronologically ordered time separated images, t₁ to t₆. A 4×6array of imaging units is marked on each time separated image. In thisexample, five (5) imaging units, IU1 to IU5, have been selected. As isclear from comparing the entire series of time separated images, t₁ tot₆, imaging unit IU1 is selected at the same position of the array (i.e.column 1, row 2) for each of the time separated images. Similarly, eachof imaging units IU2 to IU5 is selected in corresponding positions ofthe imaging unit array throughout the series of time separated images.Accordingly, values for each imaging unit at its corresponding positionthroughout the time separated images can be plotted as a function oftime to obtain a temporal curve. The plotting step is illustrated by thearrows linking the corresponding positions of imaging units IU1, IU2 andIU5. The temporal curve for each imaging unit can then be analyzed todetermine a temporal difference. For example, a mean curve for all ofthe temporal curves obtained can be calculated. The mean curve is thensubtracted from each of the temporal curves to obtain a temporaldifference curve for each selected imaging unit. It will be understood,that the marked array and the selection of imaging units IU1 to IU5, aswell as the arrows representing plotting of values for imaging unitsIU1, IU2, and IU5 as a function of time, is for illustration purposesonly, and that any number of imaging units may be selected and analyzedas desired. Furthermore, the determination of temporal differences bysubtraction of a mean curve from the temporal curves for imaging unitsis for illustration only, and other methods of analyzing the temporalcurves, for example using derivatives and/or statistical techniques arecontemplated.

A mathematical basis for blind signal separation statistical analysis ofa plurality of time separated images is now provided using timeseparated images from a CT Perfusion study as an example.

Given a series of n dynamic (i.e. time separated) images from a CTPerfusion study, the n dynamic images can be represented compactly as amatrix:{tilde over (X)}=[{tilde over (x)} ₁ ^(T) {tilde over (x)} ₂ ^(T) . . .{tilde over (x)} _(p-1) ^(T) {tilde over (x)} _(p) ^(T)]^(T)where p is the number of pixels in a CT image (i.e. p=512×512, 256×256,etc) and {tilde over (x)}_(i), i=1, . . . , p is a n×1 vector whoseelements are the pixel values in the n images or the time vs.enhancement curve from the i^(th) pixel. Note that the (i,j)^(th)element of {tilde over (X)} is the j^(th) element x_(ij) of {tilde over(x)}_(i). Each {tilde over (x)}_(i), i=1, . . . , p can also beinterpreted as repeated observations of a single time vs. enhancementcurve {tilde over (χ)} (a random process consisting of n randomvariables, χ₁, χ₂, . . . , χ_(n-1), χ_(n)).Notation is interpreted as follows:(i) a vector is represented by a tilde above a Greek or Roman alphabetcharacter with elements arranged in a column;(ii) a matrix is represented by tilde over a capitalized Greek or Romanalphabet character;(iii) the transpose of a vector or a matrix is indicated by asuperscript T on the symbol; and(iv) a random variable is represented by a Greek alphabet characterwhile an observation (samples) of a random variable is represented bythe corresponding Roman alphabet character; the correspondence betweenGreek and Roman alphabet character will be used as shown below in thefollowing discussion:χ

xψ

yζ

z(v) a random process consisting of n random variables is represented bya n×1 vector.

One approach to removing noise from the n dynamic images is to find an‘expansion’ of the p×n matrix, {tilde over (X)} and then truncate theexpansion by certain pre-determined criteria. A common expansion is thesingular value decomposition of {tilde over (X)}:{tilde over (X)}=Ũ·{tilde over (L)}·{tilde over (V)} ^(T)  (1)where Ũ is a p×p orthogonal matrix, {tilde over (V)} is a n×n orthogonalmatrix, and {tilde over (L)} is a p×n diagonal matrix with singularvalues of {tilde over (X)} as its diagonal elements. Let l₁, . . . ,l_(n) be the singular values of {tilde over (X)} (it is assumed withoutloss of generality that {tilde over (X)} is of full column rank, n andp>n). Then:{tilde over (L)}=└{tilde over (l)} ₁ . . . {tilde over (l)} _(n)┘where {tilde over (l)}_(i) i=1, . . . , n is a p×1 vector with zeroelements except for the i^(th) element which is l_(i). Usually thesingular values l₁, . . . , l_(n) are arranged in descending order andsmall singular values can be discarded as arising from noise. A simpleway to reduce noise in {tilde over (X)} is to keep only r<n singularvalues. Then:{tilde over (L)}=└{tilde over (l)} ₁ . . . {tilde over (l)} _(T){tildeover (0)}_(r+1) . . . {tilde over (0)}_(n)┘If the truncated series of singular values is used to reconstitute{tilde over (X)} according to Equation (1), then:

$\overset{\sim}{X} \approx {\begin{bmatrix}{\overset{\sim}{u}}_{1} & {\overset{\sim}{u}}_{2} & \ldots & {\overset{\sim}{u}}_{p}\end{bmatrix} \cdot \begin{bmatrix}{\overset{\sim}{1}}_{1} & \ldots & {\overset{\sim}{1}}_{r} & {\overset{\sim}{0}}_{r + 1} & \ldots & {\overset{\sim}{0}}_{n}\end{bmatrix} \cdot \begin{bmatrix}{\overset{\sim}{v}}_{1}^{T} \\\vdots \\{\overset{\sim}{v}}_{r}^{T} \\{\overset{\sim}{v}}_{r + 1}^{T} \\\vdots \\{\overset{\sim}{v}}_{n}^{T}\end{bmatrix}}$where ũ_(i) i=1, . . . , p is the i^(th) column of Ũ and {tilde over(v)}_(i) i=1, . . . , n is the i^(th) column of {tilde over (V)}.Evaluating the expansion, yields:

$\begin{matrix}\begin{matrix}{\overset{\sim}{X} \approx {\begin{bmatrix}{\overset{\sim}{u}}_{1} & {\overset{\sim}{u}}_{2} & \ldots & {\overset{\sim}{u}}_{p}\end{bmatrix} \cdot \begin{bmatrix}{\overset{\sim}{l}}_{1} & \ldots & {\overset{\sim}{1}}_{r} & {\overset{\sim}{0}}_{r + 1} & \ldots & {\overset{\sim}{0}}_{n}\end{bmatrix} \cdot \begin{bmatrix}{\overset{\sim}{v}}_{1}^{T} \\\vdots \\{\overset{\sim}{v}}_{r}^{T} \\{\overset{\sim}{v}}_{r + 1}^{T} \\\vdots \\{\overset{\sim}{v}}_{n}^{T}\end{bmatrix}}} \\{\approx {\begin{bmatrix}{\overset{\sim}{u}}_{1} & {\overset{\sim}{u}}_{2} & \ldots & {\overset{\sim}{u}}_{p}\end{bmatrix} \cdot \left\lbrack {{{\overset{\sim}{l}}_{1} \cdot {\overset{\sim}{v}}_{1}^{T}} + \ldots + {{\overset{\sim}{1}}_{r} \cdot {\overset{\sim}{v}}_{r}^{T}}} \right\rbrack}} \\{\approx {\begin{bmatrix}{\overset{\sim}{u}}_{1} & {\overset{\sim}{u}}_{2} & \ldots & {\overset{\sim}{u}}_{p}\end{bmatrix} \cdot \left( {{\begin{bmatrix}l_{1} \\\vdots \\0_{r} \\0_{r + 1} \\\vdots \\0_{p}\end{bmatrix} \cdot {\overset{\sim}{v}}_{i}^{T}} + \ldots + {\begin{bmatrix}0 \\\vdots \\l_{r} \\0_{r + 1} \\\vdots \\0_{p}\end{bmatrix} \cdot {\overset{\sim}{v}}_{r}^{T}}} \right)}} \\{\approx {\begin{bmatrix}{\overset{\sim}{u}}_{1} & {\overset{\sim}{u}}_{2} & \ldots & {\overset{\sim}{u}}_{p}\end{bmatrix} \cdot \begin{bmatrix}{l_{1} \cdot {\overset{\sim}{v}}_{1}^{T}} \\\vdots \\{l_{r} \cdot {\overset{\sim}{v}}_{r}^{T}} \\0_{r + 1} \\\vdots \\0_{p}\end{bmatrix}}} \\{= {\overset{r}{\sum\limits_{i = 1}}{l_{i} \cdot {\overset{\sim}{u}}_{i} \cdot {\overset{\sim}{v}}_{i}^{T}}}}\end{matrix} & (2)\end{matrix}$Thus, in the expansion of {tilde over (X)}, keeping r singular valuesmeans that only the first r columns of Ũ and {tilde over (V)} is needed.

The singular value decomposition (SVD) (expansion) of {tilde over (X)}as formulated here, however, has two main disadvantages in de-noisingdynamic images from CT Perfusion, namely:

1. The size of the matrix {tilde over (X)} can be as large as(256)²×40-60, making it computationally prohibitive to find the SVD of{tilde over (X)}; and

2. Although the criterion of discarding small singular values makesgeneral sense, it is difficult to justify the threshold used to keep ordiscard singular values.

An alternative method to achieve the SVD of {tilde over (X)} will now bedescribed. As stated above, the time vs. enhancement curves from eachpixel (or pixel block) {tilde over (x)}_(i), i=1, . . . , p can beregarded as repeated samples (observations) of the underlying randomprocess (time vs. enhancement curve). {tilde over (χ)}. Then the samplemean of the random process, x, can be calculated as:

$\begin{matrix}{\overset{\_}{x} = {\frac{1}{p} \cdot {\sum\limits_{i = 1}^{p}{\overset{\sim}{x}}_{i}}}} & (3)\end{matrix}$x by definition is a n×1 vector (the sample mean time vs. enhancementcurve) and the j^(th) component of x, according to Equation (3), is:

${\overset{\_}{x}}_{j} = {\frac{1}{p} \cdot {\sum\limits_{i = 1}^{p}{\overset{\sim}{x}}_{ij}}}$where {tilde over (x)}_(ij) is the j^(th) component of {tilde over(x)}_(i)The zero mean random process from is {tilde over (χ)} is {tilde over(ψ)}={tilde over (χ)}− χ, where χ is the population mean of {tilde over(χ)}. For each observation (sample) {tilde over (x)}_(i) of {tilde over(χ)}, the corresponding observation (sample) of {tilde over (ψ)} isgiven by:{tilde over (y)} _(i) ={tilde over (x)} _(i) − x.The sample mean of {tilde over (ψ)}, y, a n×1 vector is given by:

$\overset{\_}{y} = {{\frac{1}{p} \cdot {\sum\limits_{i = 1}^{p}{\overset{\sim}{y}}_{i}}} = {{\frac{1}{p} \cdot {\sum\limits_{i = 1}^{p}\left( {{\overset{\sim}{x}}_{i} - \overset{\_}{x}} \right)}} = {{\left( {\frac{1}{p} \cdot {\overset{p}{\sum\limits_{i = 1}}{\overset{\sim}{x}}_{i}}} \right) - \overset{\_}{x}} = \overset{\sim}{0}}}}$where {tilde over (0)} is a n×1 zero vector as expected since {tildeover (ψ)} is a zero mean random process.

Let {tilde over (Y)}=[{tilde over (y)}₁ ^(T) {tilde over (y)}₂ ^(T) . .. {tilde over (y)}_(p-1) ^(T) {tilde over (y)}_(p) ^(T)]^(T) be a p×nmatrix whose rows are {tilde over (y)}_(i) ^(T) i=1, . . . , p. Then:

$\begin{matrix}\begin{matrix}{{{\overset{\sim}{Y}}^{T} \cdot \overset{\sim}{Y}} = {\begin{bmatrix}{\overset{\sim}{y}}_{1} & {\overset{\sim}{y}}_{2} & \ldots & \ldots & {\overset{\sim}{y}}_{p - 1} & {\overset{\sim}{y}}_{p}\end{bmatrix} \cdot \begin{bmatrix}{\overset{\sim}{y}}_{1}^{T} \\{\overset{\sim}{y}}_{2}^{T} \\\vdots \\\vdots \\{\overset{\sim}{y}}_{p - 1}^{T} \\{\overset{\sim}{y}}_{p}^{T}\end{bmatrix}}} \\{= {\overset{p}{\sum\limits_{i = 1}}{{\overset{\sim}{y}}_{i} \cdot {\overset{\sim}{y}}_{i}^{T}}}} \\{= {\overset{p}{\sum\limits_{i = 1}}\begin{bmatrix}{\overset{\sim}{y}}_{i\; 1}^{2} & {{\overset{\sim}{y}}_{i\; 1} \cdot {\overset{\sim}{y}}_{i\; 2}} & \ldots & \ldots & {{\overset{\sim}{y}}_{i\; 1} \cdot {\overset{\sim}{y}}_{{in} - 1}} & {{\overset{\sim}{y}}_{i\; 1} \cdot {\overset{\sim}{y}}_{in}} \\{{\overset{\sim}{y}}_{i\; 2} \cdot {\overset{\sim}{y}}_{i\; 1}} & {\overset{\sim}{y}}_{i\; 2}^{2} & \ldots & \ldots & {{\overset{\sim}{y}}_{i\; 2} \cdot {\overset{\sim}{y}}_{{in} - 1}} & {{\overset{\sim}{y}}_{i\; 2} \cdot {\overset{\sim}{y}}_{in}} \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\{{\overset{\sim}{y}}_{{i\; n} - 1} \cdot {\overset{\sim}{y}}_{i\; 1}} & {{\overset{\sim}{y}}_{{in} - 1} \cdot {\overset{\sim}{y}}_{i\; 2}} & \; & \; & {\overset{\sim}{y}}_{{in} - 1}^{2} & {{\overset{\sim}{y}}_{{in} - 1} \cdot {\overset{\sim}{y}}_{in}} \\{{\overset{\sim}{y}}_{in} \cdot {\overset{\sim}{y}}_{i\; 1}} & {{\overset{\sim}{y}}_{in} \cdot {\overset{\sim}{y}}_{i\; 2}} & \; & \; & {{\overset{\sim}{y}}_{in} \cdot {\overset{\sim}{y}}_{{in} - 1}} & {\overset{\sim}{y}}_{in}^{2}\end{bmatrix}}} \\{= \begin{bmatrix}\sigma_{1}^{2} & \sigma_{12}^{2} & \ldots & \ldots & \sigma_{{in} - 1}^{2} & \sigma_{1n}^{2} \\\sigma_{21}^{2} & \sigma_{2}^{2} & \ldots & \ldots & \sigma_{{2p} - 1}^{2} & \sigma_{2n}^{2} \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\\sigma_{n - 11}^{2} & \sigma_{n - 12}^{2} & \ldots & \ldots & \sigma_{n - 1}^{2} & \sigma_{n - {1n}}^{2} \\\sigma_{n\; 1}^{2} & \sigma_{n\; 2}^{2} & \ldots & \ldots & \sigma_{{nn} - 1}^{2} & \sigma_{n}^{2}\end{bmatrix}}\end{matrix} & (4)\end{matrix}$where σ_(i) ² is the sample variance of ψ_(i) i=1, . . . , n and σ_(ij)² is the sample covariance of and ψ_(i) and ψ_(j) i,j=1, . . . , n.

Or,

$\overset{\sim}{S} = {\frac{1}{p - 1} \cdot {\overset{\sim}{Y}}^{T} \cdot \overset{\sim}{Y}}$is the sample covariance matrix of {tilde over (ψ)} and is of dimensionn×n. In CT Perfusion n is typically 40-60.A linear combination of ψ_(i), i=1, . . . , n, ζ, can be defined as:ζ=ã ^(T)·{tilde over (ψ)}where ã is a given n×1 vector of coefficients. ζ is a random variable(not process) with observations z_(i) i=1, . . . , p given by:z _(i) =ã ^(T) ·{tilde over (y)} _(i)by definition of ζ. The sample mean of ζ, z, is given by:

$\overset{\_}{z} = {{\frac{1}{p} \cdot {\sum\limits_{i = 1}^{p}{\overset{\sim}{z}}_{i}}} = {{\frac{1}{p} \cdot {\sum\limits_{i = 1}^{p}{{\overset{\sim}{a}}^{T} \cdot {\overset{\sim}{y}}_{i}}}} = {{{\overset{\sim}{a}}^{T} \cdot \frac{1}{p} \cdot {\sum\limits_{i = 1}^{p}{\overset{\sim}{y}}_{i}}} = 0}}}$Thus, ζ is also a zero mean random variable. The sample variance of ζ isgiven by:

$\begin{matrix}{{{sample}\mspace{14mu}{{var}(\zeta)}} = {\frac{1}{p - 1} \cdot {\sum\limits_{i = 1}^{p}\left( z_{i} \right)^{2}}}} \\{= {\frac{1}{p - 1} \cdot {\sum\limits_{i = 1}^{p}{\left( {{\overset{\sim}{a}}^{T} \cdot {\overset{\sim}{y}}_{i}} \right) \cdot \left( {{\overset{\sim}{a}}^{T} \cdot {\overset{\sim}{y}}_{i}} \right)^{T}}}}} \\{= {\frac{1}{p - 1} \cdot {\sum\limits_{i = 1}^{p}{{\overset{\sim}{a}}^{T} \cdot {\overset{\sim}{y}}_{i} \cdot {\overset{\sim}{y}}_{i}^{T} \cdot \overset{\sim}{a}}}}} \\{= {\frac{1}{p - 1} \cdot {\overset{\sim}{a}}^{T} \cdot \left( {\sum\limits_{i = 1}^{p}{\cdot {\overset{\sim}{y}}_{i} \cdot {\overset{\sim}{y}}_{i}^{T}}} \right) \cdot \overset{\sim}{a}}} \\{= {{\overset{\sim}{a}}^{T} \cdot \overset{\sim}{S} \cdot \overset{\sim}{a}}}\end{matrix}$

Determination of the vector of coefficients, ã, to maximize the samplevariance of ζ, ã^(T)·{tilde over (S)}·ã a will now be investigated.Since ã^(T)·{tilde over (S)}·ã can be made as large as required byscaling ã with a constant, the optimization can only be done by imposinga normalization condition, such as ã^(T)·ã=1. To maximize ã^(T)·{tildeover (S)}·ã subject to ã^(T)·ã=1, the standard approach is to use thetechnique of Lagrange multipliers. That is it is sought to maximize:ã ^(T) ·{tilde over (S)}·ã−λ·(ã ^(T) ·ã−1)where λ is the Lagrange multiplier. Differentiation with respect to ãgives:{tilde over (S)}·ã−λ·ã={tilde over (0)}Thus, λ is the eigenvalue of {tilde over (S)} and ã is the correspondingeigenvector andsample var(ζ)=ã ^(T) ·{tilde over (S)}·ã=ã ^(T) ·λ·ã=λ, the eigenvalueof {tilde over (S)}That is to maximize ã^(T)·{tilde over (S)}·ã subject to ã^(T)·ã=1, ã isthe eigenvector of {tilde over (S)} that has the highest eigenvalue.

Let λ₁ be the largest eigenvalue with the corresponding eigenvector ã₁,then ζ₁=ã₁ ^(T)·{tilde over (ψ)} is called the first sample principalcomponent of {tilde over (ψ)}. It is a random variable with observationsdefined by z_(i1)=ã₁ ^(T)·{tilde over (y)}_(i)·z_(i1) is also called thescore of the i^(th) observation on the first principal component. Thesecond principal component, ζ₂=ã₂ ^(T)·{tilde over (ψ)}, similarly,maximizes the sample variance of ζ₂, ã₂ ^(T)·{tilde over (S)}·ã₂,subject to being uncorrelated to ζ₁ and the normalization conditional ã₂^(T)·ã₂=1. The observations corresponding to ζ₂ is z_(i2)=α₂ ^(T)·{tildeover (y)}_(i). The sample covariance of ζ₁ and ζ₂ is:

$\begin{matrix}{{{sample}\mspace{14mu}{{cov}\left( {\zeta_{1},\zeta_{2}} \right)}} = {\frac{1}{p - 1} \cdot {\sum\limits_{i = 1}^{p}{z_{i\; 1} \cdot z_{i\; 2}}}}} \\{= {\frac{1}{p - 1} \cdot {\sum\limits_{i = 1}^{p}{\left( {{\overset{\sim}{a}}_{1}^{T} \cdot {\overset{\sim}{y}}_{i}} \right) \cdot \left( {{\overset{\sim}{a}}_{2}^{T} \cdot {\overset{\sim}{y}}_{i}} \right)^{T}}}}} \\{= {\frac{1}{p - 1} \cdot {\sum\limits_{i = 1}^{p}{{\overset{\sim}{a}}_{1}^{T} \cdot {\overset{\sim}{y}}_{i} \cdot {\overset{\sim}{y}}_{i}^{T} \cdot {\overset{\sim}{a}}_{2}}}}} \\{= {\frac{1}{p - 1} \cdot {\overset{\sim}{a}}_{1}^{T} \cdot \left( {\sum\limits_{i = 1}^{p}{\cdot {\overset{\sim}{y}}_{i} \cdot {\overset{\sim}{y}}_{i}^{T}}} \right) \cdot {\overset{\sim}{a}}_{2}}} \\{= {{\overset{\sim}{a}}_{1}^{T} \cdot \overset{\sim}{S} \cdot {\overset{\sim}{a}}_{2}}} \\{= {\frac{1}{p - 1} \cdot {\sum\limits_{i = 1}^{p}{\left( {{\overset{\sim}{a}}_{2}^{T} \cdot {\overset{\sim}{y}}_{i}} \right) \cdot \left( {{\overset{\sim}{a}}_{1}^{T} \cdot {\overset{\sim}{y}}_{i}} \right)^{T}}}}} \\{= {{\overset{\sim}{a}}_{2}^{T} \cdot \overset{\sim}{S} \cdot {\overset{\sim}{a}}_{1}}}\end{matrix}$The requirement that ζ₂ is uncorrelated to ζ₁ is the same as requiring:ã ₂ ^(T) ·{tilde over (S)}·ã ₁ =ã ₂ ^(T)·λ₁ ·ã ₁=λ₁ ·[ã ₂ ^(T) ·ã ₁]=0

ã ₂ ^(T) ·{tilde over (S)}·ã ₁=0 and ã ₂ ^(T) ·ã ₁=0ã ₁ ^(T) ·{tilde over (S)}·ã ₂=({tilde over (S)}·ã ₁)^(T) ·ã ₂=λ₁ ·[ã ₁^(T) ·ã ₂]=0

ã ₁ ^(T) ·{tilde over (S)}·ã ₂=0 and ã ₁ ^(T) ·ã ₂=0

Thus, any of the conditions:ã ₂ ^(T) ·{tilde over (S)}·ã ₁=0ã ₂ ^(T) ·ã ₁=0ã ₁ ^(T) ·{tilde over (S)}·ã ₂=0ã ₁ ^(T) ·ã ₂=0can be used to specify that ζ₂ is uncorrelated to ζ₁.

To maximize the sample variance of ζ₂, ã₂ ^(T)·{tilde over (S)}·ã₂,subject to the conditions that ζ₂ is uncorrelated to ζ₁ and ã₂ ^(T)·ã₂=1is the same as to maximize:ã ₂ ^(T) ·{tilde over (S)}·ã ₂−λ₂·(ã ₂ ^(T) ·ã ₂−1)−φã ₂ ^(T) ·ã ₁where the condition ã₂ ^(T)·ã₁=0 is used to specify that ζ₂ isuncorrelated to ζ₁ and λ₂ and φ are Lagrange multipliers.Differentiation with respect to ã₂ gives:{tilde over (S)}·ã ₂−λ₂ ·ã ₂ −φ·ã ₁=0Multiplying through by ã₁ ^(T) gives:ã ₁ ^(T) ·{tilde over (S)}·ã ₂−λ₂ ·ã ₁ ^(T) ·ã ₂ −φ·ã ₁ ^(T) ·ã ₁=0which, since the first two terms are zero and ã₁ ^(T)·ã₁=1, reduces toφ=0. Therefore,{tilde over (S)}·ã ₂−λ₂ ·ã ₂=0so λ₂ is once more an eigenvalue of {tilde over (S)}, and ã₂ thecorresponding eigenvector.

Again, ã₂ ^(T)·{tilde over (S)}·ã₂=λ₂ or the sample variance of ζ₂ isλ₂. Assuming that {tilde over (S)} does not have repeat eigenvalues, λ₂cannot be equal to λ₁. If it did, it follows that ã₂=ã₁, violating theconstraint that ã₂ ^(T)·ã₁=0. Hence λ₂ is the second largest eigenvalueand ã₂ the corresponding eigenvector.

It can be shown using the same technique that for the third, fourth, . .. , n^(th) sample principal components, the vectors of coefficients ã₃,ã₄, . . . , ã_(n) are the eigenvectors of {tilde over (S)} correspondingto λ₃, λ₄, . . . , λ_(n), the third, fourth largest, . . . , and thesmallest eigenvalue. Furthermore, the sample variances of the principalcomponents are also given by the eigenvalues of {tilde over (S)}.

Define {tilde over (Z)} to be the p×n matrix of the scores ofobservations on principal components. In this case, the i^(th) rowconsists of the scores of the i^(th) observation of {tilde over (ψ)}(i.e. {tilde over (y)}_(i)) on all principal components ζ₁, ζ₂, . . . ,ζ_(n-1), ζ_(n), that is:

$\begin{matrix}\begin{matrix}{\overset{\sim}{Z} = \begin{bmatrix}{{\overset{\sim}{a}}_{1}^{T} \cdot {\overset{\sim}{y}}_{1}} & \ldots & {{\overset{\sim}{a}}_{j}^{T} \cdot {\overset{\sim}{y}}_{1}} & \ldots & {{\overset{\sim}{a}}_{n}^{T} \cdot {\overset{\sim}{y}}_{1}} \\\vdots & \ddots & \vdots & \ddots & \vdots \\{{\overset{\sim}{a}}_{1}^{T} \cdot {\overset{\sim}{y}}_{i}} & \ldots & {{\overset{\sim}{a}}_{j}^{T} \cdot {\overset{\sim}{y}}_{i}} & \ldots & {{\overset{\sim}{a}}_{n}^{T} \cdot {\overset{\sim}{y}}_{i}} \\\vdots & \ddots & \vdots & \ddots & \vdots \\{{\overset{\sim}{a}}_{1}^{T} \cdot {\overset{\sim}{y}}_{p}} & \ldots & {{\overset{\sim}{a}}_{j}^{T} \cdot {\overset{\sim}{y}}_{p}} & \ldots & {{\overset{\sim}{a}}_{n}^{T} \cdot {\overset{\sim}{y}}_{p}}\end{bmatrix}} \\{= \begin{bmatrix}{{\overset{\sim}{a}}_{1}^{T} \cdot {\overset{\sim}{y}}_{1}} & \ldots & {{\overset{\sim}{a}}_{1}^{T} \cdot {\overset{\sim}{y}}_{i}} & \ldots & {{\overset{\sim}{a}}_{1}^{T} \cdot {\overset{\sim}{y}}_{p}} \\\vdots & \ddots & \vdots & \ddots & \vdots \\{{\overset{\sim}{a}}_{j}^{T} \cdot {\overset{\sim}{y}}_{1}} & \ldots & {{\overset{\sim}{a}}_{j}^{T} \cdot {\overset{\sim}{y}}_{i}} & \ldots & {{\overset{\sim}{a}}_{j}^{T} \cdot {\overset{\sim}{y}}_{p}} \\\vdots & \ddots & \vdots & \ddots & \vdots \\{{\overset{\sim}{a}}_{n}^{T} \cdot {\overset{\sim}{y}}_{1}} & \ldots & {{\overset{\sim}{a}}_{n}^{T} \cdot {\overset{\sim}{y}}_{i}} & \ldots & {{\overset{\sim}{a}}_{n}^{T} \cdot {\overset{\sim}{y}}_{p}}\end{bmatrix}^{T}} \\{= \left( {\begin{bmatrix}{\overset{\sim}{a}}_{1}^{T} \\\; \\{\overset{\sim}{a}}_{j}^{T} \\\; \\{\overset{\sim}{a}}_{n}^{T}\end{bmatrix} \cdot \left\lbrack {{\overset{\sim}{y}}_{1}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{y}}_{i}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{y}}_{p}} \right\rbrack} \right)^{T}} \\{= {\left\lbrack {{\overset{\sim}{y}}_{1}^{T}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{y}}_{i}^{T}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{y}}_{p_{p}}^{T}} \right\rbrack^{T} \cdot \left\lbrack {{\overset{\sim}{a}}_{1}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{a}}_{j}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{a}}_{n}} \right\rbrack}} \\{= {\overset{\sim}{Y} \cdot \overset{\sim}{A}}}\end{matrix} & \left( {4a} \right)\end{matrix}$where Ã is a n×n matrix whose columns are the eigenvectors of {tildeover (S)} and is orthogonal.

The first principal component, ζ₁, maximized is expressed by:

${{\overset{\sim}{a}}_{1}^{T} \cdot \overset{\sim}{S} \cdot {\overset{\sim}{a}}_{1}} = {{\frac{1}{p - 1} \cdot {\overset{\sim}{a}}_{1}^{T} \cdot \left( {\sum\limits_{i = 1}^{p}{\cdot {\overset{\sim}{y}}_{i} \cdot {\overset{\sim}{y}}_{i}^{T}}} \right) \cdot {\overset{\sim}{a}}_{1}} = {\frac{1}{p - 1} \cdot {\overset{\sim}{a}}_{1}^{T} \cdot \left( {\sum\limits_{i = 1}^{p}{\left( {{\overset{\sim}{x}}_{i} - \overset{\_}{x}} \right) \cdot \left( {{\overset{\sim}{x}}_{i} - \overset{\_}{x}} \right)^{T}}} \right) \cdot {\overset{\sim}{a}}_{1}}}$({tilde over (x)}_(i)− x) is the deviation of the time vs. enhancementcurve (TEC) from the i^(th) pixel (pixel block) from the mean of TECsfrom all pixels (pixel blocks). Thus, the eigenvectors: ã₁, ã₂, ã₃, ã₄,. . . , ã_(n) represent the first, second, third, fourth, . . . , andleast dominant time vs. enhancement behaviour in the set of TECs fromall pixels. The loadings of the eigenvectors ã₁, ã₂, ã₃, ã₄, . . . ,ã_(n) in the TEC from the i^(th) pixel (pixel block) are:ã ₁ ^(T) ·{tilde over (y)} _(i) ,ã ₂ ^(T) ·{tilde over (y)} _(i) ,ã ₃^(T) ·{tilde over (y)} _(i) ,ã ₄ ^(T) ·{tilde over (y)} _(i) , . . . ,ã_(n) ^(T) ·{tilde over (y)} _(i)and is the i^(th) row of the matrix {tilde over (Z)} while the loadingsof the eigenvector ã_(i) in the TECs from all pixels (pixel blocks) are:ã _(i) ^(T) ·{tilde over (y)} ₁ ,ã _(i) ^(T) ·{tilde over (y)} ₂ ,ã _(i)^(T) ·{tilde over (y)} ₃ ,ã _(i) ^(T) ·{tilde over (y)} ₄ , . . . , ã_(i) ^(T) ·{tilde over (y)} _(p)and is the i^(th) column of the matrix {tilde over (Z)}. The i^(th)column of {tilde over (Z)}, therefore, gives the loading map of theeigenvector ã_(i) which corresponds to the i^(th) largest eigenvalueλ_(i).

To de-noise the dynamic images from a CT Perfusion study, the reverse ofthe calculation of the loading maps of eigenvectors is required.Reconstituting a smooth version of the original dynamic images, {tildeover (Y)}({tilde over (X)}), from a truncated series of eigenvectors ã₁,ã₂, ã₃, ã₄, . . . , ã_(r) is necessary.

As in Equation (1), the p×n matrix {tilde over (Y)} can be expanded bysingular value decomposition as:{tilde over (Y)}=Ũ·{tilde over (L)}·{tilde over (V)} ^(T)where Ũ is a p×p orthogonal matrix, {tilde over (V)} is a n×n orthogonalmatrix, and {tilde over (L)} is a p×n diagonal matrix with singularvalues of {tilde over (Y)} as its diagonal elements. Let l₁, . . . ,l_(n) be the singular values of {tilde over (Y)} (it is assumed withoutloss of generality that {tilde over (Y)} is of full column rank, n andp>n). Using the singular value decomposition:

$\begin{matrix}{\overset{\sim}{S} = {\frac{1}{p - 1} \cdot {\overset{\sim}{Y}}^{T} \cdot \overset{\sim}{Y}}} \\{= {\frac{1}{p - 1} \cdot \left( {\overset{\sim}{U} \cdot \overset{\sim}{L} \cdot {\overset{\sim}{V}}^{T}} \right)^{T} \cdot \overset{\sim}{U} \cdot \overset{\sim}{L} \cdot {\overset{\sim}{V}}^{T}}} \\{= {\frac{1}{p - 1} \cdot \overset{\sim}{V} \cdot {\overset{\sim}{L}}^{T} \cdot {\overset{\sim}{U}}^{T} \cdot \overset{\sim}{U} \cdot \overset{\sim}{L} \cdot {\overset{\sim}{V}}^{T}}} \\{= {\frac{1}{p - 1} \cdot \overset{\sim}{V} \cdot {\overset{\sim}{L}}^{T} \cdot \overset{\sim}{L} \cdot {\overset{\sim}{V}}^{T}}}\end{matrix}$Multiplying both sides of the equation by {tilde over (V)} yields:

$\begin{matrix}{{\overset{\sim}{S} \cdot \overset{\sim}{V}} = {\frac{1}{p - 1} \cdot \overset{\sim}{V} \cdot {\overset{\sim}{L}}^{T} \cdot \overset{\sim}{L}}} & (5)\end{matrix}${tilde over (L)}=└l₁ . . . {tilde over (l)}_(i) . . . {tilde over(l)}_(n)┘ where {tilde over (l)}_(i) i=1, . . . , n is a p×1 vectorwhose elements are all zero except the i^(th) element which is equal tol_(i), that is {tilde over (l)}_(i) ^(T)=[0 . . . l_(i) . . . 0] and

${{\overset{\sim}{1}}_{i}^{T} \cdot {\overset{\sim}{1}}_{j}} = \left\{ {{\begin{matrix}\begin{matrix}1_{i}^{2} & {i = j}\end{matrix} \\\begin{matrix}0 & {i \neq j}\end{matrix}\end{matrix}{Then}},\mspace{14mu}{{{\overset{\sim}{L}}^{T} \cdot \overset{\sim}{L}} = {{\begin{bmatrix}{\overset{\sim}{1}}_{1}^{T} \\\vdots \\{\overset{\sim}{1}}_{i}^{T} \\\vdots \\{\overset{\sim}{1}}_{n}^{T}\end{bmatrix} \cdot \left\lbrack {{\overset{\sim}{1}}_{1}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{1}}_{i}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{1}}_{n}} \right\rbrack} = \begin{bmatrix}{\overset{\sim}{m}}_{1}^{T} \\\vdots \\{\overset{\sim}{m}}_{i}^{T} \\\vdots \\{\overset{\sim}{m}}_{n}^{T}\end{bmatrix}}}} \right.$where {tilde over (m)}_(i) ^(T)=[0 . . . l_(i) ² . . . 0].Let {tilde over (V)}=[{tilde over (v)}₁ . . . {tilde over (v)}_(i) . . .{tilde over (v)}_(n)] where {tilde over (v)}_(i) i=1, . . . , n is a n×1vector.Then,

${\overset{\sim}{V} \cdot {\overset{\sim}{L}}^{T} \cdot \overset{\sim}{L}} = {{\left\lbrack {{\overset{\sim}{v}}_{1}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{v}}_{i}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{v}}_{n}} \right\rbrack \cdot \begin{bmatrix}{\overset{\sim}{m}}_{1}^{T} \\\vdots \\{\overset{\sim}{m}}_{i}^{T} \\\vdots \\{\overset{\sim}{m}}_{n}^{T}\end{bmatrix}} = {{\sum\limits_{i = 1}^{n}{{\overset{\sim}{v}}_{i} \cdot {\overset{\sim}{m}}_{i}^{T}}} = {{\sum\limits_{i = 1}^{n}{{\overset{\sim}{v}}_{i} \cdot \left\lbrack {0\mspace{14mu}\ldots\mspace{14mu} 1_{i}^{2}\mspace{14mu}\ldots\mspace{14mu} 0} \right\rbrack}} = {{\sum\limits_{i = 1}^{n}\left\lbrack {0\mspace{14mu}\ldots\mspace{14mu}{1_{i}^{2} \cdot {\overset{\sim}{v}}_{i}}\mspace{14mu}\ldots\mspace{14mu} 0} \right\rbrack} = \left\lbrack {{1_{1}^{2} \cdot {\overset{\sim}{v}}_{1}}\mspace{14mu}\ldots\mspace{14mu}{1_{i}^{2} \cdot {\overset{\sim}{v}}_{i}}\mspace{14mu}\ldots\mspace{14mu}{1_{n}^{2} \cdot {\overset{\sim}{v}}_{n}}} \right\rbrack}}}}$Substituting in Equation (5) yields:

${{\overset{\sim}{S} \cdot \overset{\sim}{V}} = {{\frac{1}{p - 1} \cdot \overset{\sim}{V} \cdot {\overset{\sim}{L}}^{T} \cdot \overset{\sim}{L}} = {\left. \left\lbrack {\frac{1_{1}^{2} \cdot {\overset{\sim}{v}}_{1}}{p - 1}\mspace{14mu}\ldots\mspace{14mu}\frac{1_{i}^{2} \cdot {\overset{\sim}{v}}_{i}}{p - 1}\mspace{14mu}\ldots\mspace{14mu}\frac{1_{n}^{2} \cdot {\overset{\sim}{v}}_{n}}{p - 1}} \right\rbrack\Rightarrow\left\lbrack {{\overset{\sim}{S} \cdot {\overset{\sim}{v}}_{1}}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{S} \cdot {\overset{\sim}{v}}_{i}}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{S} \cdot {\overset{\sim}{v}}_{n}}} \right\rbrack \right. = {\left. \left\lbrack {\frac{1_{1}^{2} \cdot {\overset{\sim}{v}}_{1}}{p - 1}\mspace{14mu}\ldots\mspace{14mu}\frac{1_{i}^{2} \cdot {\overset{\sim}{v}}_{i}}{p - 1}\mspace{14mu}\ldots\mspace{14mu}\frac{1_{n}^{2} \cdot {\overset{\sim}{v}}_{n}}{p - 1}} \right\rbrack\Rightarrow{\overset{\sim}{S} \cdot {\overset{\sim}{v}}_{1}} \right. = {{{\frac{1_{i}^{2}}{p - 1} \cdot {\overset{\sim}{v}}_{i}}\mspace{31mu} i} = 1}}}}},\ldots\mspace{14mu},n$From the discussion on eigenvectors and eigenvalues:

${\overset{\sim}{v}}_{i} = {\left. {\overset{\sim}{a}}_{i}\Rightarrow\overset{\sim}{V} \right. = \overset{\sim}{A}}$$\frac{1_{i}^{2}}{p - 1} = \lambda_{i}$or the singular values of {tilde over (Y)}(l_(i)) is equal to the squareroot of (p−1) times the eigenvalue value of {tilde over (S)}(λ_(i)),where

${\overset{\sim}{S} = {\frac{1}{p - 1} \cdot {\overset{\sim}{Y}}^{T} \cdot \overset{\sim}{Y}}},$i.e. l_(i)=√{square root over ((p−1)·λ_(i))}Thus,{tilde over (Y)} ^(T) ·{tilde over (Y)}·{tilde over (v)} _(i) =l _(i) ²·{tilde over (v)} _(i) or {tilde over (Y)} ^(T) ·{tilde over (Y)}·ã _(i)=l _(i) ² ·ã _(i) i=1, . . . , n  (5a)Next, consider

$\begin{matrix}\begin{matrix}{{\overset{\sim}{S}}^{*} = {\frac{1}{p - 1} \cdot \overset{\sim}{Y} \cdot {\overset{\sim}{Y}}^{T}}} \\{= {\frac{1}{p - 1} \cdot \overset{\sim}{U} \cdot \overset{\sim}{L} \cdot {\overset{\sim}{V}}^{T} \cdot \left( {\overset{\sim}{U} \cdot \overset{\sim}{L} \cdot {\overset{\sim}{V}}^{T}} \right)^{T}}} \\{= {\frac{1}{p - 1}{\overset{\sim}{U} \cdot \overset{\sim}{L} \cdot {\overset{\sim}{V}}^{T} \cdot \overset{\sim}{V} \cdot \overset{\sim}{L^{T}} \cdot {\overset{\sim}{U}}^{T}}}} \\{= {\frac{1}{p - 1} \cdot \overset{\sim}{U} \cdot L \cdot {\overset{\sim}{L}}^{T} \cdot {\overset{\sim}{U}}^{T}}}\end{matrix} & (6)\end{matrix}$(note that {tilde over (S)}* is not the transpose of {tilde over (S)})Multiplying both sides of the equation by Ũ yields:

$\mspace{20mu}{{{\overset{\sim}{S}}^{*} \cdot \overset{\sim}{U}} = {\frac{1}{p - 1} \cdot \overset{\sim}{U} \cdot \overset{\sim}{L} \cdot {\overset{\sim}{L}}^{T}}}$$\mspace{20mu}{{\overset{\sim}{L} \cdot {\overset{\sim}{L}}^{T}} = {{\left\lbrack {{\overset{\sim}{1}}_{1}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{1}}_{i}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{1}}_{n}} \right\rbrack \cdot \begin{bmatrix}{\overset{\sim}{1}}_{1}^{T} \\\vdots \\{\overset{\sim}{1}}_{i}^{T} \\\vdots \\{\overset{\sim}{1}}_{n}^{T}\end{bmatrix}} = {\sum\limits_{i = 1}^{n}{{\overset{\sim}{1}}_{1} \cdot {\overset{\sim}{1}}_{i}^{T}}}}}$${{\overset{\sim}{1}}_{1} \cdot {\overset{\sim}{1}}_{i}^{T}} = {{\underset{\underset{p \times 1}{︸}}{\begin{bmatrix}0_{1} \\\vdots \\1_{i} \\\vdots \\0_{p}\end{bmatrix}} \cdot \underset{\underset{1 \times p}{︸}}{\left\lbrack {0_{1}\mspace{14mu}\ldots\mspace{14mu} 1_{i}\mspace{14mu}\ldots\mspace{14mu} 0_{p}} \right\rbrack}} = \begin{bmatrix}0 & \ldots & 0_{i} & \ldots & 0_{1n} & \ldots & 0_{1p} \\\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\0_{i\; 1} & \ldots & 1_{i}^{2} & \ldots & 0_{i\; n} & \ldots & 0_{ip} \\\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\0_{n\; 1} & \ldots & 0_{ni} & \ldots & 0_{nn} & \ldots & 0_{np} \\\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\0_{p\; 1} & \ldots & 0_{pi} & \ldots & 0_{pn} & \ldots & 0_{pp}\end{bmatrix}}$$\mspace{20mu}{{{\overset{\sim}{1}}_{i} \cdot {\overset{\sim}{1}}_{j}^{T}} = {{\underset{\underset{p \times 1}{︸}}{\begin{bmatrix}0_{1} \\\vdots \\1_{i} \\\vdots \\0_{j} \\\vdots \\0_{p}\end{bmatrix}} \cdot \underset{\underset{1 \times p}{︸}}{\left\lbrack {0_{1}\mspace{14mu}\ldots\mspace{14mu} 0_{i}\mspace{14mu}\ldots\mspace{14mu} 1_{j}\mspace{14mu}\ldots\mspace{14mu} 0_{p}} \right\rbrack}} = {\overset{\sim}{0}}_{p \times p}}}$$\mspace{20mu}{{Thus},\mspace{20mu}{{\overset{\sim}{L} \cdot {\overset{\sim}{L}}^{T}} = {\begin{bmatrix}1_{1}^{2} & \ldots & 0_{i} & \ldots & 0_{1n} & \ldots & 0_{1p} \\\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\0_{i\; 1} & \ldots & 1_{i}^{2} & \ldots & 0_{i\; n} & \ldots & 0_{ip} \\\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\0_{n\; 1} & \ldots & 0_{ni} & \ldots & 1_{n}^{2} & \ldots & 0_{np} \\\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\0_{p\; 1} & \ldots & 0_{pi} & \ldots & 0_{pn} & \ldots & 0_{pp}\end{bmatrix} = \begin{bmatrix}{\overset{\sim}{m}}_{1}^{*T} \\\vdots \\{\overset{\sim}{m}}_{i}^{*T} \\\vdots \\{\overset{\sim}{m}}_{n}^{*T} \\\vdots \\0\end{bmatrix}}}}$Let {tilde over (m)}_(i)*=[0 . . . l_(i) ² . . . 0_(p)]^(T) be a p×1vector as defined, then

${\overset{\sim}{L} \cdot {\overset{\sim}{L}}^{T}} = \begin{bmatrix}{\overset{\sim}{m}}_{1}^{*T} \\\vdots \\{\overset{\sim}{m}}_{i}^{*T} \\\vdots \\{\overset{\sim}{m}}_{n}^{*T} \\0_{n + 1} \\\vdots \\0_{p}\end{bmatrix}$Let Ũ=[ũ₁ . . . ũ_(i) . . . ũ_(p)] where ũ_(i) i=1, . . . , p is a p×1vector.Then,

$\begin{matrix}{{\overset{\sim}{U} \cdot \overset{\sim}{L} \cdot {\overset{\sim}{L}}^{T}} = {\left\lbrack {{\overset{\sim}{u}}_{1}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{u}}_{i}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{u}}_{n}\mspace{14mu}{\overset{\sim}{u}}_{n + 1}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{u}}_{p}} \right\rbrack \cdot \begin{bmatrix}{\overset{\sim}{m}}_{1}^{*T} \\\vdots \\{\overset{\sim}{m}}_{i}^{*T} \\\vdots \\{\overset{\sim}{m}}_{n}^{*T} \\0_{n + 1} \\\vdots \\0_{p}\end{bmatrix}}} \\{= {\sum\limits_{i = 1}^{n}{{{\overset{\sim}{u}}_{i} \cdot {\overset{\sim}{m}}_{i}^{T}}1}}} \\{= {{\sum\limits_{i = 1}^{n}{{\overset{\sim}{u}}_{i} \cdot \left\lbrack {0\mspace{14mu}\ldots\mspace{14mu} 1_{i}^{2}\mspace{14mu}\ldots\mspace{14mu} 0_{p}} \right\rbrack}} + {\sum\limits_{i = {n + 1}}^{p}{{\overset{\sim}{u}}_{i} \cdot \left\lbrack {0\mspace{14mu}\ldots\mspace{14mu} 0_{i}\mspace{14mu}\ldots\mspace{14mu} 0_{p}} \right\rbrack}}}} \\{= {\sum\limits_{i = 1}^{n}\left\lbrack {0\mspace{14mu}\ldots\mspace{14mu}{1_{i}^{2} \cdot {\overset{\sim}{u}}_{i}}\mspace{14mu}\ldots\mspace{14mu} 0_{p}} \right\rbrack}} \\{= \left\lbrack {{1_{1}^{2} \cdot {\overset{\sim}{u}}_{i}}\mspace{14mu}\ldots\mspace{14mu}{1_{i}^{2} \cdot {\overset{\sim}{u}}_{i}}\mspace{14mu}\ldots\mspace{14mu}{1_{n}^{2} \cdot {\overset{\sim}{u}}_{n}}\mspace{14mu} 0_{n + 1}\mspace{14mu}\ldots\mspace{14mu} 0_{p}} \right\rbrack}\end{matrix}$Substituting in Equation (6) yields:

${{\overset{\sim}{S}}^{*} \cdot \overset{\sim}{U}} = {{\frac{1}{p - 1} \cdot \overset{\sim}{U} \cdot \overset{\sim}{L} \cdot {\overset{\sim}{L}}^{T}} = {\left. \left\lbrack {\frac{1_{1}^{2} \cdot {\overset{\sim}{u}}_{1}}{p - 1}\mspace{14mu}\ldots\mspace{14mu}\frac{1_{i}^{2} \cdot {\overset{\sim}{u}}_{i}}{p - 1}\mspace{14mu}\ldots\mspace{14mu}\frac{1_{n}^{2} \cdot {\overset{\sim}{u}}_{n}}{p - 1}\mspace{14mu} 0_{n + 1}\mspace{14mu}\ldots\mspace{14mu} 0_{p}} \right\rbrack\Longrightarrow\left\lbrack {{{\overset{\sim}{S}}^{*} \cdot {\overset{\sim}{u}}_{i}}\mspace{14mu}\ldots\mspace{14mu}{{\overset{\sim}{S}}^{*} \cdot {\overset{\sim}{u}}_{i}}\mspace{14mu}\ldots\mspace{14mu}{{\overset{\sim}{S}}^{*} \cdot {\overset{\sim}{u}}_{n}}\mspace{14mu} 0_{n + 1}\mspace{14mu}\ldots\mspace{14mu} 0_{p}} \right\rbrack \right. = {\quad\left\lbrack {{\frac{1_{1}^{2} \cdot {\overset{\sim}{u}}_{1}}{p - 1}\mspace{14mu}\ldots\mspace{14mu}\frac{1_{i}^{2} \cdot {\overset{\sim}{u}}_{i}}{p - 1}\mspace{14mu}\ldots\mspace{14mu}\frac{1_{n}^{2} \cdot {\overset{\sim}{u}}_{n}}{p - 1}\mspace{14mu}{\left. 0_{n + 1}\mspace{20mu}\Longrightarrow{\overset{\sim}{S}}^{*} \right. \cdot {\overset{\sim}{u}}_{i}}} = \left\{ \begin{matrix}{\frac{1_{i}^{2}}{p - 1} \cdot {\overset{\sim}{u}}_{i}} & {i \leq n} \\{0 \cdot {\overset{\sim}{u}}_{i}} & {{n + 1} \leq i \leq p}\end{matrix} \right.} \right.}}}$Thus, ũ_(i) i=1, . . . , n is the eigenvector of {tilde over (S)}*corresponding to the i^(th) eigenvalue

$\frac{1_{i}^{2}}{p - 1},$which happens also to be the i^(th) eigenvalue of {tilde over (S)},while all ũ_(i) i=n+1, . . . , p are eigenvectors of {tilde over (S)}*with zero eigenvalue. Also{tilde over (Y)}·{tilde over (Y)} ^(T) ·ũ _(i) =l _(i) ² ·ũ _(i) i=1, .. . , n  (6a)As stated above, Equation (5a) is:{tilde over (Y)} ^(T) ·{tilde over (Y)}·{tilde over (v)} _(i) =l _(i) ²·{tilde over (v)} _(i)Multiplying both sides by {tilde over (Y)} yields:{tilde over (Y)}·{tilde over (Y)} ^(T) ·{tilde over (Y)}·{tilde over(v)} _(i) =l _(i) ² ·{tilde over (Y)}·{tilde over (v)} _(i)Comparing with Equation (6a):{tilde over (Y)}·{tilde over (v)} _(i) =c·ũ _(i) where c is aconstant  (7a)Similarly multiplying both sides of Equation (6a) by {tilde over(Y)}^(T) yields:{tilde over (Y)} ^(T) ·{tilde over (Y)}·{tilde over (Y)} ^(T) ·ũ _(i) =l_(i) ² ·{tilde over (Y)} ^(T) ·ũ _(i)Comparing with Equation (5a):{tilde over (Y)} ^(T) ·ũ _(i) =c*·{tilde over (v)} _(i) where c is aconstant  (7b)From Equation (7b):

${\overset{\sim}{v}}_{i} = {\frac{1}{c^{*}} \cdot {\overset{\sim}{Y}}^{T} \cdot {\overset{\sim}{u}}_{i}}$Substituting into Equation (7a) yields:{tilde over (Y)}·{tilde over (Y)} ^(T) ·ũ _(i) =c*·c·ũ _(i)Substituting into Equation (6a) yields:l _(i) ² ·ũ _(i) =c*·c·ũ _(i)Therefore, c*·c=l_(i) ²Without loss of generality, c* is set to:c*=c=l _(i)=√{square root over ((p−1)·λ_(i))} i=1, . . . , n and n<pThus, for the singular value decomposition of {tilde over (Y)}:{tilde over (Y)}=Ũ·{tilde over (L)}·{tilde over (V)} ^(T)where Ũ is a p×p orthogonal matrix:Ũ=[ũ ₁ . . . ũ _(i) . . . ũ _(p)] and ũ _(i) i=1, . . . , p is a p×1vector;{tilde over (V)} is a n×n orthogonal matrix:{tilde over (V)}=[{tilde over (v)} ₁ . . . {tilde over (v)} _(i) . . .{tilde over (v)} _(n)] and {tilde over (v)} _(i) i=1, . . . , n is a n×1vector;and {tilde over (L)} is a p×n diagonal matrix with singular values of{tilde over (Y)}, l₁, . . . , l_(n), as its diagonal elements; it isalso assumed without loss of generality that {tilde over (Y)} is of fullcolumn rank, n and p>n.It is shown that:

{tilde over (v)}_(i)=ã_(i) where ã_(i) is the i^(th) eigenvector of

$\overset{\sim}{S} = {\frac{1}{p - 1} \cdot {\overset{\sim}{Y}}^{T} \cdot \overset{\sim}{Y}}$corresponding to the eigenvalue

$\frac{1_{i}^{2}}{\left( {p - 1} \right)};$and from Equation (7c)

${\overset{\sim}{u}}_{i} = {{\frac{1}{c}{\overset{\sim}{Y} \cdot {\overset{\sim}{v}}_{i}}} = {\frac{1}{\sqrt{\left( {p - 1} \right) \cdot \lambda_{i}}} \cdot \overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{i}}}$As discussed above, to de-noise {tilde over (Y)} the SVD expansion canbe truncated after the r^(th) singular value:

$\begin{matrix}{{\overset{\sim}{Y} \approx {\sum\limits_{i = 1}^{r}{1_{i} \cdot {\overset{\sim}{u}}_{i} \cdot {\overset{\sim}{v}}_{i}^{T}}}} = {{\sum\limits_{i = 1}^{r}{1_{i} \cdot {\overset{\sim}{u}}_{i} \cdot {\overset{\sim}{a}}_{i}^{T}}} = {{\sum\limits_{i = 1}^{r}{1_{i} \cdot \frac{1}{\sqrt{\left( {p - 1} \right) \cdot \lambda_{i}}} \cdot \overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{i} \cdot {\overset{\sim}{a}}_{i}^{T}}} = {{\sum\limits_{i = 1}^{r}{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{i} \cdot \left. {\overset{\sim}{a}}_{i}^{T}\Longrightarrow\overset{\sim}{Y} \right.}} \approx {\overset{\sim}{Y} \cdot {\sum\limits_{i = 1}^{r}{{\overset{\sim}{a}}_{i} \cdot {\overset{\sim}{a}}_{i}^{T}}}}}}}} & (8)\end{matrix}$

Equation (8) also affords a very simple geometric interpretation asexplained in the following.

$\begin{matrix}{{\overset{\sim}{Y} \approx {\overset{\sim}{Y} \cdot {\sum\limits_{i = 1}^{r}{{\overset{\sim}{a}}_{i} \cdot {\overset{\sim}{a}}_{i}^{T}}}}} = {{\sum\limits_{i = 1}^{r}{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{i} \cdot {\overset{\sim}{a}}_{i}^{T}}} = {{\left( {{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{1}}\mspace{14mu}{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{2}}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{r - 1}}\mspace{14mu}{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{r}}} \right) \cdot \begin{pmatrix}{\overset{\sim}{a}}_{1}^{T} \\{\overset{\sim}{a}}_{2}^{T} \\\vdots \\{\overset{\sim}{a}}_{r - 1}^{T} \\{\overset{\sim}{a}}_{r}^{T}\end{pmatrix}} = {\overset{\sim}{Y} \cdot {\overset{\sim}{A}}_{r} \cdot {\overset{\sim}{A}}_{r}^{T}}}}} & (9)\end{matrix}$where Ã_(r) is a n×r matrix whose i^(th) column is ã_(i) and fromEquation (4a):{tilde over (Z)} _(r) ={tilde over (Y)}·Ã _(r)and the rows and columns of {tilde over (Z)}_(r) have the sameinterpretation as those of {tilde over (Z)}.Rewrite Equation (9):

$\begin{matrix}{{\overset{\sim}{Y} \approx {\overset{\sim}{Y} \cdot {\sum\limits_{i = 1}^{r}{{\overset{\sim}{a}}_{i} \cdot {\overset{\sim}{a}}_{i}^{T}}}}} = {{\sum\limits_{i = 1}^{r}{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{i} \cdot {\overset{\sim}{a}}_{i}^{T}}} = {{\left( {{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{1}}\mspace{14mu}{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{2}}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{r - 1}}\mspace{14mu}{\overset{\sim}{Y} \cdot {\overset{\sim}{a}}_{r}}} \right) \cdot \begin{pmatrix}{\overset{\sim}{a}}_{1}^{T} \\{\overset{\sim}{a}}_{2}^{T} \\\vdots \\{\overset{\sim}{a}}_{r - 1}^{T} \\{\overset{\sim}{a}}_{r}^{T}\end{pmatrix}} = {{\left( {{\begin{pmatrix}{\overset{\sim}{y}}_{1}^{T} \\{\overset{\sim}{y}}_{2}^{T} \\\vdots \\{\overset{\sim}{y}}_{p - 1}^{T} \\{\overset{\sim}{y}}_{p}^{T}\end{pmatrix} \cdot {\overset{\sim}{a}}_{1}}\mspace{14mu}{\begin{pmatrix}{\overset{\sim}{y}}_{1}^{T} \\{\overset{\sim}{y}}_{2}^{T} \\\vdots \\{\overset{\sim}{y}}_{p - 1}^{T} \\{\overset{\sim}{y}}_{p}^{T}\end{pmatrix} \cdot {\overset{\sim}{a}}_{2}}\mspace{14mu}\ldots\mspace{14mu}{\begin{pmatrix}{\overset{\sim}{y}}_{1}^{T} \\{\overset{\sim}{y}}_{2}^{T} \\\vdots \\{\overset{\sim}{y}}_{p - 1}^{T} \\{\overset{\sim}{y}}_{p}^{T}\end{pmatrix} \cdot {\overset{\sim}{a}}_{r - 1}}\mspace{14mu}{\begin{pmatrix}{\overset{\sim}{y}}_{1}^{T} \\{\overset{\sim}{y}}_{2}^{T} \\\vdots \\{\overset{\sim}{y}}_{p - 1}^{T} \\{\overset{\sim}{y}}_{p}^{T}\end{pmatrix} \cdot {\overset{\sim}{a}}_{r}}} \right) \cdot \begin{pmatrix}{\overset{\sim}{a}}_{1}^{T} \\{\overset{\sim}{a}}_{2}^{T} \\\vdots \\{\overset{\sim}{a}}_{r - 1}^{T} \\{\overset{\sim}{a}}_{r}^{T}\end{pmatrix}} = {{\begin{pmatrix}{{\overset{\sim}{y}}_{1}^{T} \cdot {\overset{\sim}{a}}_{1}} & {{\overset{\sim}{y}}_{1}^{T} \cdot {\overset{\sim}{a}}_{2}} & \ldots & {{\overset{\sim}{y}}_{1}^{T} \cdot {\overset{\sim}{a}}_{r - 1}} & {{\overset{\sim}{y}}_{1}^{T} \cdot {\overset{\sim}{a}}_{r}} \\{{\overset{\sim}{y}}_{2}^{T} \cdot {\overset{\sim}{a}}_{1}} & {{\overset{\sim}{y}}_{2}^{T} \cdot {\overset{\sim}{a}}_{2}} & \ldots & {{\overset{\sim}{2}}_{2}^{T} \cdot {\overset{\sim}{a}}_{r - 1}} & {{\overset{\sim}{2}}_{2}^{T} \cdot {\overset{\sim}{a}}_{r}} \\\vdots & \vdots & \ldots & \vdots & \vdots \\{{\overset{\sim}{y}}_{p - 1}^{T} \cdot {\overset{\sim}{a}}_{1}} & {{\overset{\sim}{y}}_{p - 1}^{T} \cdot {\overset{\sim}{a}}_{2}} & \ldots & {{\overset{\sim}{y}}_{p - 1}^{T} \cdot {\overset{\sim}{a}}_{r - 1}} & {{\overset{\sim}{y}}_{p - 1}^{T} \cdot {\overset{\sim}{a}}_{r}} \\{{\overset{\sim}{y}}_{p}^{T} \cdot {\overset{\sim}{a}}_{1}} & {{\overset{\sim}{y}}_{p}^{T} \cdot {\overset{\sim}{a}}_{2}} & \ldots & {{\overset{\sim}{y}}_{p}^{T} \cdot {\overset{\sim}{a}}_{r - 1}} & {{\overset{\sim}{y}}_{p}^{T} \cdot {\overset{\sim}{a}}_{r}}\end{pmatrix} \cdot \begin{pmatrix}{\overset{\sim}{a}}_{1}^{T} \\{\overset{\sim}{a}}_{2}^{T} \\\vdots \\{\overset{\sim}{a}}_{r - 1}^{T} \\{\overset{\sim}{a}}_{r}^{T}\end{pmatrix}} = \begin{pmatrix}{\sum\limits_{i = 1}^{r}{\left( {{\overset{\sim}{y}}_{1}^{T} \cdot {\overset{\sim}{a}}_{i}} \right) \cdot {\overset{\sim}{a}}_{i}^{T}}} \\{\sum\limits_{i = 1}^{r}{\left( {{\overset{\sim}{y}}_{2}^{T} \cdot {\overset{\sim}{a}}_{i}} \right) \cdot {\overset{\sim}{a}}_{i}^{T}}} \\\vdots \\{\sum\limits_{i = 1}^{r}{\left( {{\overset{\sim}{y}}_{p - 1}^{T} \cdot {\overset{\sim}{a}}_{i}} \right) \cdot {\overset{\sim}{a}}_{i}^{T}}} \\{\sum\limits_{i = 1}^{r}{\left( {{\overset{\sim}{y}}_{p}^{T} \cdot {\overset{\sim}{a}}_{i}} \right) \cdot {\overset{\sim}{a}}_{i}^{T}}}\end{pmatrix}}}}}} & (10)\end{matrix}$

Equation (10) suggests that the TEC for the i^(th) pixel (pixel block)is a weighted summation of r eigenvectors with the weights given by theloading of the i^(th) pixel TEC on the eigenvectors.

The eigenvectors ã₁, ã₂, . . . , ã_(r−1), ã_(r) can be regarded asforming an orthonormal basis in n dimensional space. The projection of{tilde over (y)}_(i) on ã_(i) is {tilde over (y)}_(i) ^(T)·ã_(i) and{tilde over (y)}_(i) can be reconstituted as

$\sum\limits_{j = 1}^{r}{\left( {{\overset{\sim}{y}}_{i}^{T} \cdot {\overset{\sim}{a}}_{j}} \right) \cdot {{\overset{\sim}{a}}_{j}^{T}.}}$

The algorithm to de-noise {tilde over (X)}=[{tilde over (x)}₁ ^(T) . . .{tilde over (x)}_(i) ^(T) . . . {tilde over (x)}_(p) ^(T)] is asfollows:

1. Calculate the p×n zero mean matrix, {tilde over (Y)}=[{tilde over(y)}₁ ^(T) . . . {tilde over (y)}_(i) ^(T) . . . {tilde over (y)}_(p)^(T)]^(T), corresponding to {tilde over (X)}, where

${\overset{\sim}{y}}_{i} = {{{\overset{\sim}{x}}_{i} - {\overset{\_}{x}\mspace{14mu}{and}\mspace{14mu}\overset{\_}{x}}} = {\frac{1}{p} \cdot {\sum\limits_{i = 1}^{p}{\overset{\sim}{x}}_{i}}}}$2. Calculate the n×n sample covariance matrix, {tilde over (S)}:

$\overset{\sim}{S} = {\frac{1}{p - 1} \cdot {\overset{\sim}{Y}}^{T} \cdot \overset{\sim}{Y}}$3. Calculate the n pairs of eigenvalue and eigenvectors,(λ_(i) ,ã _(i))i=1, . . . , n4. Calculate the n×n smoothing matrix by keeping only the largest reigenvalues:

$\sum\limits_{i = 1}^{r}{{\overset{\sim}{a}}_{i} \cdot {\overset{\sim}{a}}_{i}^{T}}$5. Smooth {tilde over (Y)}=[{tilde over (y)}₁ . . . {tilde over (y)}_(i). . . {tilde over (y)}_(p)] by calculating the product:

${\overset{\sim}{Y}}^{*} = {\left\lbrack {{\overset{\sim}{y}}_{1}^{*}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{y}}_{i}^{*}\mspace{14mu}\ldots\mspace{14mu}{\overset{\sim}{y}}_{p}^{*}} \right\rbrack \approx {\overset{\sim}{Y} \cdot {\sum\limits_{i = 1}^{r}{{\overset{\sim}{a}}_{i} \cdot {\overset{\sim}{a}}_{i}^{T}}}}}$6. Calculate a smoothed version of {tilde over (X)}* as follows:{tilde over (X)}*=[{tilde over (y)} ₁ *+ x . . . {tilde over (y)} _(i)*+ x . . . {tilde over (y)} _(p) *+ x].

An example of the method described herein will now be demonstrated withregard to experimental results. The principal component analysis (PCA)described herein is tested for the ability to achieve dose reduction ofCT scans.

The PCA method will now be summarized based on a CT Perfusion study thatgenerates a plurality of time separated images that capture the passageof contrast through the brain. Corresponding to each pixel in the seriesof images, a curve is generated that represents the contrastconcentration in that pixel as a function of time. A difference curve,equal to the difference of the curve from the mean of all pixel curvesis also generated. A covariance matrix of the co-variations of all pairsof difference curves is calculated and represents all possiblevariations about the mean curve. The principal component analysis methodanalyzes the covariance matrix to find common temporal patterns amongthe curves by finding linear combinations of the difference curves thatrepresent most of the variations about the mean curve. It can be shownthat such linear combinations, or principal components, are theeigenvectors of the covariance matrix and that the eigenvector with thelargest eigenvalue is the most prominent temporal feature among all thedifference curves, the eigenvector with the second largest eigenvalue isthe second most prominent temporal feature and so on, while eigenvectorswith small eigenvalues are from noise in the CT images. A corollary isthat each difference curve can be represented as a sum of all theprincipal components and by discarding components with smalleigenvalues, noise can be effectively suppressed. To denoise the timeseries of images, the mean curve is added to all the difference curvesreconstituted by using the first few principal components with largeeigenvalues.

While currently available statistical filters remove noise in projectiondata before image reconstruction, the principal component analysismethod removes noise from the reconstructed images. Experimental resultsdescribed below show that performing PCA leads to dose reduction in CTimaging.

The PCA method is tested versus the traditional filtered back-projection(FBP) method in radiation dose reduction for CT Perfusion applications.The dose saving of PCA versus FBP on measuring brain perfusion using CTPerfusion was tested in a healthy pig, which was scanned four times withthe x-ray tube current changed from 190 mA to 100 mA; 75 mA and 50 mA,respectively, while keeping other scanning parameters the same.

FIG. 3 shows brain blood flow maps of the pig calculated using the J-Wmodel with CT Perfusion images obtained using (A) 190 mA and FBPreconstruction; (B) 100 mA, FBP reconstruction and PCA; (C) 70 mA, FBPreconstruction and PCA; and, (D) 50 mA, FBP reconstruction and PCA. FIG.3(E) shows regions of interest (3 outlined circles) used for calculationof Figure of Merit of the blood flow maps shown in FIG. 4. The regionsare superimposed on a CT image which shows a coronal section through thebrain of the pig. The dark areas are the skull and jaw bones.

FIG. 4 shows that in terms of standard deviation/mean of blood flow(Figure of Merit) in the three regions shown in FIG. 3(E): 70 mA withPCA was better than 190 mA with FBP, while 50 mA PCA was worse than 190mA with FBP (p<0.05). The image quality of the blood flow maps in FIG. 3also reflects this progression. Therefore, the PCA method is able toreduce radiation dose required to obtain images in CT Perfusion studiesand in maintaining the quality of derived blood flow maps compared tothe use of FBP alone. More specifically, the PCA method achievedapproximately three fold reduction in radiation dose compared to the useof FBP alone. The effect of the PCA method on dose reduction inproduction of BBB-PS maps was not demonstrated because the blood-brainbarrier remained intact. However, since blood flow, blood volume andBBB-PS are calculated together in the J-W model, the dose reductionshown using blood flow maps is expected to apply to BBB-PS maps.

The principal component analysis methodology described herein has beenimplemented in a computer executable software application using C++. Thetime required to process ninety-five (95) 512×512 CT brain images isless than two (2) minutes on a general purpose computing device such asfor example a personal computer (PC).

The method described herein can be implemented on hardware such assystem 400 shown in FIG. 5. An input 510 receives a plurality of timeseparated images, which can be previously stored, received directly froman imaging device, or communicated over a wired or wirelesscommunication link from a remote location. A general purpose computingdevice 505 receives the time separated images from the input andexecutes a software application thereby to process the time separatedimages as described above. The general purpose computing device 505interfaces with the user through a display 520, a keyboard 515 and/or amouse or other pointing device 525. Results, for example the constructedprocessed images, can be presented on display 520 and/or output to anysuitable output device 530, for example, a printer, a storage device, ora communication device for communicating the results to a remotelocation.

The software application described herein may run as a stand-aloneapplication or may be incorporated into other available applications toprovide enhanced functionality to those applications. The softwareapplication may include program modules including routines, programs,object components, data structures etc. and may be embodied as computerreadable program code stored on a computer readable medium. The computerreadable medium is any data storage device that can store data, whichcan thereafter be read by a computer system. Examples of computerreadable media include for example read-only memory, random-accessmemory, CD-ROMs, magnetic tape and optical data storage devices. Thecomputer readable program code can also be distributed over a networkincluding coupled computer systems so that the computer readable programcode is stored and executed in a distributed fashion.

The above-described embodiments are intended to be examples andalterations and modifications may be effected thereto, by those of skillin the art, without departing from the scope of the invention which isdefined by the claims appended hereto.

What is claimed is:
 1. A method of medical image processing comprising:capturing a plurality of time separated images of a target region ofinterest of a patient, wherein the plurality of time separated imagescapture a time-varying behavior of said target region of interest;selecting a plurality of imaging units from each of the plurality ofcaptured time separated images, wherein the selected plurality ofimaging units are at the same location within each of the plurality ofcaptured time separated images; representing each of the selectedplurality of imaging units by a value; for each imaging unit location,using the corresponding values to determine a temporal difference forthe selected imaging units; selecting temporal differences above athreshold limit; and constructing a noise reduced image based on theselected temporal differences.
 2. The method of claim 1, furthercomprising ranking each measured temporal difference prior to saidselecting.
 3. The method of claim 1, wherein each imaging unit is one ofa pixel and a voxel.
 4. The method of claim 1, wherein each imaging unitis one of a plurality of pixels and a plurality of voxels.
 5. The methodof claim 1, wherein the temporal difference is based on contrastenhancement.
 6. The method of claim 5, wherein the temporal differenceis a difference of a contrast agent.
 7. The method of claim 6, whereincontrast enhancement is based on the concentration of the contrastagent.
 8. The method of claim 6, wherein contrast enhancement is basedon the activation of the contrast agent.
 9. The method of claim 1,further comprising characterizing the measured temporal differences witha blind signal separation statistical technique.
 10. The method of claim1, wherein determining the temporal difference comprises calculating acovariance matrix of temporal differences for the imaging units.
 11. Themethod of claim 10, further comprising analyzing the covariance matrixwith a blind signal separation statistical technique to determine commontemporal patterns among the temporal differences.
 12. The method ofclaim 11, further comprising analyzing the eigenvectors of thecovariance matrix.
 13. The method of claim 12, wherein eigenvalues ofthe eigenvector represent the degree of temporal difference.
 14. Themethod of claim 13, wherein the selecting of temporal differencescomprises selecting eigenvalues above the threshold limit.
 15. Themethod of claim 1, wherein the threshold limit is adjustable.
 16. Themethod of claim 11, wherein the blind signal separation statisticaltechnique is selected from the group consisting of principal componentanalysis, independent component analysis, blind source deconvolution,Karhunen-Loeve transform and Hotelling transform.
 17. The method ofclaim 1, wherein the plurality of time separated images are capturedusing a dynamic medical imaging technique.
 18. The method of claim 17,wherein the dynamic medical imaging technique comprises use of acontrast agent.
 19. The method of claim 17, wherein the time separatedimages are captured using a dynamic medical imaging technique thatemploys a protocol having a reduced radiation dose.
 20. The method ofclaim 17, wherein the time separated images are captured using a dynamicmedical imaging technique that employs a protocol having a reducedamount of contrast agent.
 21. The method of claim 17, wherein thedynamic medical imaging technique is magnetic resonance imaging (MRI),functional MRI, positron emission tomography, nuclear medicine orcomputed tomography.
 22. The method of claim 17, wherein the dynamicmedical imaging technique is computed tomography (CT).
 23. The method ofclaim 22, wherein the time separated images are captured using a CTPerfusion protocol having a reduced radiation dose.
 24. The method ofclaim 22, wherein the time separated images are captured using a CTPerfusion protocol having a reduced amount of contrast agent.
 25. Asystem for medical image processing comprising: an interface configuredto receive a plurality of captured time separated images of a targetregion of interest of a patient, wherein the plurality of captured timeseparated images show a time-varying behavior of said target region ofinterest; and a processing structure communicating with said interfaceand configured to process said plurality of captured time separatedimages, wherein during processing of said captured time separatedimages, said processing structure is configured to select a plurality ofimaging units from each of the plurality of captured time separatedimages, wherein the selected imaging units are at the same locations ineach of said plurality of captured time separated images, represent eachof the selected imaging units by a value, for each imaging unitlocation, use the values representing the selected imaging units in eachof the time separated images to determine a temporal difference for theselected imaging units, select temporal differences above a thresholdlimit, and construct a noised reduced image based on the selectedtemporal differences.
 26. A non-transitory computer readable mediumembodying a computer program for medical image processing, the computerprogram comprising computer program code, which when executed by one ormore processors, causes an apparatus at least to: capture a plurality oftime separated images of a target region of interest of a patient,wherein the plurality of captured time separated images image atime-varying behavior of said target region of interest; select aplurality of imaging units from each of the plurality of captured timeseparated images, wherein the selected imaging units are at the samelocations in each of said captured images; represent each of theselected imaging units by a value; for each imaging unit location, usethe values representing the selected imaging units in each of the timeseparated images to determine a temporal difference for the selectedimaging units; select temporal differences above a threshold limit; andconstruct a noise reduced image based on the selected temporaldifferences.